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Abstract. In this paper, we present an image separation method for 
separating images into point- and curvehke parts by employing a com- 
bined dictionary consisting of wavelets and compactly supported shear- 
lets utilizing the fact that they sparsely represent point and curvilinear 
singularities, respectively. Our methodology is based on the very recently 
introduced mathematical theory of geometric separation, which shows 
that highly precise separation of the morphologically distinct features of 
points and curves can be achieved by i'^ minimization. Finally, we present 
some experimental results showing the effectiveness of our algorithm, in 
particular, the ability to accurately separate points from curves even if 
the curvature is relatively large due to the excellent localization property 
of compactly supported shearlets. 

Keywords: Geometric separation, £i minimization, sparse approxima- 
tion, shearlets, wavelets 



1 Introduction 

The task of separating an image into its morphologically different contents has 
recently drawn a lot of attention in the research community due to its signif- 
icance for applications. In neurobiological imaging, it would, for instance, be 
desirable to separate 'spines' (pointlike objects) from 'dendrites' (curvelike ob- 
jects) in order to analyze them independently aiming to detect characteristics 
of Alzheimer disease. Also, in astronomical imaging, astronomers would often 
like to separate stars from filaments for further analysis, hence again separating 
point- from curvelike structures. Successful methodologies for efficiently and ac- 
curately solving this task can in fact be applied to a much broader range of areas 
in science and technology including medical imaging, surveillance, and speech 
processing. 

Although the problem of separating morphologically distinct features seems 
to be intractable - the problem is underdetermined, since there is only one known 
data (the image) and two or more unknowns - there has been extensive studies 
on this topic. The book by Meyer [18] initiated the area of image decomposition, 
in particular, the utilization of variational methods. Some years later, Starck, 
Elad, and Donoho suggested a different approach in [19 coined 'Morphological 
Component Analysis', which proclaims that such a separation task might be 
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possible provided that we have prior information about the type of features to 
be extracted and provided that the morphological difference between those is 
strong enough. For the separation of point- and curvelike features, it was in fact 
recently even theoretically proven in [3] that ii minimization solves this task 
with arbitrarily high precision exploring a combined dictionary of wavelets and 
curvelets. Wavelets provide optimally sparse expansions for point like structures, 
and curvelets provide optimally sparse expansions for curvelike structures. Thus 
ii minimization applied to the expansion coefficients of the original image into 
this combined dictionary forces the pointlike structures into the wavelet part and 
the curvelike structures into the curvelet part, thereby automatically separating 
the image. An associated algorithmic approach using wavelets and curvelets has 
been implemented in MCALatEl. 

Recently, a novel directional representation system - so-called shearlets - 
has emerged which provides a unified treatment of continuum models as well 
as digital models, allowing, for instance, a precise resolution of wavefront sets, 
optimally sparse representations of cartoon-like images, and associated fast de- 
composition algorithms; see the survey paper [11]. Shear let systems are systems 
generated by one single generator with parabolic scaling, shearing, and transla- 
tion operators applied to it, in the same way wavelet systems are dyadic scalings 
and translations of a single function, but including a directionality character- 
istic owing to the additional shearing operation (and the anisotropic scaling). 
The shearing operation in fact provides a more favorable treatment of direc- 
tions, thereby ensuring a unified treatment of the continuum and digital realm 
as opposed to curvelets which are rotation-based in the continuum realm, see 

Thus, it is natural to ask whether also a combined dictionary of wavelets 
and shearlets might be utilizable for separating point- and curvelike features, 
the advantage presumably being a faster scheme, a more precise separation, 
and a direct applicability of theoretical results achieved for the continuum do- 
main. And, in fact, the theoretical results from [3 based on a model situation 
were shown to also hold for a combined dictionary of wavelets and shearlets [2] . 
Moreover, numerical results give evidence to the superior behavior of shear let- 
based decomposition algorithms when compared to curvelet-based algorithms; 
see [11] for a comparison of ShearLatS with 

CurveLatEI. 

In this paper, we will present a novel approach to the separation of point- 
and curvelike features exploiting a combined dictionary of wavelets and shearlets 
as well as utilizing block relaxation in a particular way. Numerical results give 
evidence that indeed the previously anticipated advantages hold true, i.e., that 
this approach is superior to separation algorithms using wavelets and curvelets 
such as MCALab in various ways, in particular, our algorithm is faster and pro- 
vides a more precise separation, in particular, if the curvature of the curvilinear 



^ MCALab (Version 120) is available from |http : // j starck . free . f r/ j starck/Home . html [ 
^ ShearLab (Version 1.0) is available from http : //www . shear lab . org! 
^ CurveLab (Version 2.1.2) is available from http : / / www . curvelet . org] 
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part is large. In the spirit of reproducible research [5,, our algorithm is included 
in the freely available ShearLab toolbox. 

This paper is organized as follows. In Section [2l we introduce the multiscale 
system of shearlets, and Section [3] reviews the mathematical theory of geometric 
separation of point- and curvelike features. Our novel algorithmic approach is 
presented in Section |4] with numerical results discussed in Section [5l 

2 Shearlets 

In most multivariate problems, important features of the considered data are 
concentrated on lower dimensional manifolds. For example, in image processing 
an edge is an ID curve that follows a path of rapid change in image inten- 
sity. Recently, the novel directional representation system of shearlets |14|8j has 
emerged to provide efficient tools for analyzing the intrinsic geometrical features 
of a signal using anisotropic and directional window functions. In this approach, 
directionality is achieved by applying integer powers of a shear matrix, and those 
operations preserve the structure of the integer lattice which is crucial for dig- 
ital implementations. In fact, this key idea leads to a unified treatment of the 
continuum as well as digital realm, while still providing optimally sparse approx- 
imations of anisotropic features. As already mentioned before, shearlet systems 
are generated by parabolic scaling, shearing, and translation operators applied to 
one single generator. Let us now be more precise and formally introduce shearlet 
systems in 2D. 

We first start with some definitions for later use. For j > and /c G Z, let 

= ( 2?/^) ' = (^^Q 2^,) , and Sk = 

We can now define so-called cone-adapted discrete shearlet systems, where the 
term 'cone-adapted' originates from the fact that these systems tile the frequency 
domain in a cone-like fashion. For this, let c be a positive constant, which will 
later control the sampling density. For 0,?/^,?/^ G L^(R^), the cone- adapted dis- 
crete shearlet system SH{(j), ijj^ilj^c) is then defined by 

SH{(^, V^, V^; c) = ^{^; c) U ^{t. c) U #(^; c), 

where 

^((/); c) = {(/){' -cm) :me 1?}, 

^(V^;c) = {V^,,fc,^ = 25X^fcA2. • -cm) :j >0, 1^1 < [2^'/2^ , m G Z^}, 

#(^;c) = {4fc,^ = 25^'^(5ji2. • -cm) :j >0, 1^1 < \2^l%m^l?}. 

In [9^, a comprehensive theory of compactly supported shearlet frames is 
provided, i.e., systems with excellent spatial localization. It should also be men- 
tioned that in [12 a large class of compactly supported shearlet frames were 
shown to provide optimally sparse approximations of images governed by curvi- 
linear structures, in particular, so-called cartoon-like images as defined in [1]. 
This fact will be explored in the sequel. 
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3 Mathematical Theory of Geometric Separation 

In [19], a novel image separation method - Morphological Component Analy- 
sis (MCA) - based on sparse representations of images was introduced. In this 
approach, it is assumed that each image is the linear combination of several 
components that are morphologically distinct - for instance, points, curves, and 
textures. The success of this method relies on the assumption that each of the 
components is sparsely represented in a specific representation system. The key 
idea is then the following: Provided that such representation systems are iden- 
tified, the usage of a pursuit algorithm searching for the sparsest representation 
of the image with respect to the dictionary combining all those specific repre- 
sentation systems will lead to the desired separation. 

Various experimental results in [T9] show the effectiveness of this method for 
image separation however without any accompanying mathematical justification. 
Recently, the first author of this paper and Donoho developed a mathematical 
framework in [3] within which the notion of successful separation can be made 
definitionally precise and can be mathematically proven in case of separating 
point- from curvelike features, which they coined Geometric Separation. One key 
ingredient of their analysis is the consideration of clustered sparsity properties 
measured by so-called cluster coherence. In this section, we briefiy review this 
theoretical approach to the Geometric Separation Problem, which will serve as 
the foundation for our algorithm. 



3.1 Model Situation 

As a mathematical model for a composition of point- and curvelike structures, 
we consider the following two components: As a 'point-like' object, we consider 
the function V which is smooth except for point singularities and is defined by 

V = f2\x-x.\-'/\ 

As a 'curve-like' object, we consider the distribution C with singularity along a 
closed curve r : [0, 1] ^M? defined by 




Then our model situation is the sum of both, i.e., 

f = r + c. (1) 

The Geometric Separation Problem now consists of recovering V and C from the 
observed signal /. 
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3.2 Chosen Dictionary 

As we indicated before, it is now crucial to choose two representations systems 
each of which sparsely represents one of the morphologically different compo- 
nents in the Geometric Separation Problem. Our sparse approximation result, 
described in the previous section, suggests that curvilinear singularities can be 
sparsely represented by shearlets. On the other hand, it is well known that 
wavelets can provide optimally sparse approximations of functions which are 
smooth apart from point singularities. Hence, we choose the overcomplete sys- 
tem within which we will expand the signal / as a composition of the following 
two systems: 

— Orthonormal Separable Meyer Wavelets: Band-limited wavelets which form 
an orthonormal basis of isotropic generating elements. 

— Bandlimited Shearlets: A directional and anisotropic tight frame generated 
by a band-limited shearlet generator i/j defined in Section [2l 

3.3 Subband Filtering 

Since the scaling subbands of shearlets and wavelets are similar we can define 
a family of filters {Fj)j which allows to decompose a function / into pieces fj 
with different scales j depending on those subbands. The piece fj associated to 
subband j arises from filtering / using Fj by 

= F, * /, 

resulting in a function whose Fourier transform fj is supported on the scaling 
subband of scale j of the wavelet as well as the shearlet frame. The filters are 
defined in such way, that the original function can be reconstructed from the 
sequence {fj)j using 

3 

We can now exploit these tools to attack the Geometric Separation Problem 
scale-by-scale. For this, we filter the model problem ([T]) to derive the sequence 
of filtered images 

fj = Vj + Cj for all scales j. 

3.4 ^1 Minimization Problem 

Let now and ^2 be an orthonormal basis of band-limited wavelets and a tight 
frame of band-limited shearlets, respectively. Then, for each scale j, we consider 
the following optimization problem: 



{W„S,) = B.rgmm^^^sMW,h + \\1>lSih subject to f^ = W,+S,. (2) 
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Notice that ^JWj and are the wavelet and shear let coefficients of the 

signals Wj and Sj^ respectively. Notice that our objective is not on searching for 
the sparsest expansion in a wavelet-shearlet dictionary, but on separation. Thus 
we can avoid an extensive, presumably numerically instable search by minimizing 
specific coefficients, namely the analysis in contrast to the synthesis coefficients, 
for each possible separation fj = Wj + Sj . 

We wish to further remark, that here the £i norm is placed on the analysis 
rather than the synthesis coefficients to avoid numerical instabilities due to the 
redundancy of the shear let frame. 



3.5 Theoretical Result 

The theoretical result of the precision of separation of fj via (|2]) proved in |3] 
and [2] can now be stated in the following way: 

Theorem 1 ([3j and |2]). Let Wj and Sj be solutions to the optimization 
problem (0) for each scale j . Then we have 



\\rj-Wj\\2 + \\C,-Sjh 

11^.112 + IIQIb 



0, j 



This result shows that the components Vj and Cj are recovered with asymp- 
totically arbitrarily high precision at very fine scales. The energy in the point like 
component is completely captured by the wavelet coefficients, and the curvelike 
component is completely contained in the shearlet coefficients. Thus, the theory 
evidences that the Geometric Separation Problem can be satisfactorily solved 
by using a combined dictionary of wavelets and shearlets and an appropriate £i 
minimization problem. 



3.6 Extensions 

Our numerical scheme for image separation, which we will present in detail in 
the next section, will use a shift invariant wavelet tight frame and a compactly 
supported shearlet frame as opposed to orthonormal Meyer wavelets and band- 
limited shearlets required for Theorem [H Hence this deserves some comments. 
Firstly, Theorem [T] is based on an abstract separation estimate which holds for 
any pair of frames, provided certain relative sparsity and cluster coherence condi- 
tions with respect to the components of the data to be separated are satisfied (cf. 
[3]). Secondly, using the recently introduced concept of sparsity equivalence (see 
|2|10j ). results requiring sparsity and coherence conditions can be transferred 
from one system (set of systems) to another by 'merely' considering particular 
decay conditions of the cross-Grammian matrix (matrices). We strongly believe 
that this framework allows a similar result as Theorem [1] for the pair of a shift 
invariant wavelet tight frame and a compactly supported shearlet frame. Since 
the focus of this paper is however on the introduction of the numerical scheme, 
such a highly technical, theoretical analysis is beyond the scope of this paper 
and will be treated in a subsequent work. 
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4 Our Algorithmic Approach to the Geometric 
Separation Problem 

In this section, we present our algorithmic approach to the Geometric Separa- 
tion Problem of separating point- from curvelike features by using a combined 
dictionary of wavelets and shearlets. The ingredients of the algorithm will be 
detailed below. 

4.1 General Scheme 

In practice, the observed signal / is often contaminated by noise which requires 
an adaption of the optimization problem ([2]). As proposed in numerous publica- 
tions, one typically considers a modified optimization problem - so-called Basis 
Pursuit Denoising (BPDN) - which can be obtained by relaxing the constraint 
in in order to deal with noisy observed signals (see [6l). For each scale j, the 
optimization problem then takes the form: 



In this new form, the additional content in the image - the noise - charac- 
terized by the property that it can not be represented sparsely by either one of 
the two representation systems will be allocated to the residual fj — Wj — Sj. 
Hence, performing this minimization, we not only separate point- and curvelike 
objects, which were modeled by Vj and Cj in Subsection 13. ![ but also succeed in 
removing an additive noise component as a by-product. Of course, solving the 
optimization problem (|3]) for all relevant scales j is computationally expensive. 

4.2 Preprocessing 

To avoid high complexity, we observe that the frequency distribution of point- 
and curvelike components is highly concentrated on high frequencies. Hence it 
would be essentially sufficient for achieving accurate separation to solve (jSj) for 
only sufficiently large scales j, as also evidenced by Theorem [H This idea leads 
to a simplification of the problem (j3j) by modifying the observed signal / as 
follows: We first consider bandpass filters Fq, . . . ,Fz,, where (i^j)j=o,...,L is the 
family of bandpass filters defined in Subsection 13.31 up to scale L, and Fq is a 
lowpass filter. Thus the observed signal / satisfies 



For each scale j, we now carefully choose a non-uniform weight wj > satisfying, 
in particular, Wj < Wj'^ if j < j' . These weights are then utilized for a weighted 
reconstruction of / resulting in a newly constructed signal / by computing 



L 



f = T.^J^fj with/,=F,*/. 



L 




(4) 
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In this way, the two morphological components, namely points and curves, can 
be enhanced by suppressing the low frequencies. 

While emphasizing that certainly other weights can be applied in our scheme, 
for the numerical tests presented in this paper we chose L = 3, i.e., 4 subbands, 
and weights wq = O^wi = 0.1, = 0.7, and ws = 0.7. This coincides with the 
intuition that a strong weight should be assigned to a band on which the power 
spectrum of the underlying morphological contents is highly concentrated. We 
do not claim that this is necessarily the optimal choice, and a comprehensive 
mathematical optimality analysis is beyond our reach at this moment. To our 
mind, the numerical results though justify this choice. 

4.3 Solver for the £i Minimization Problem 

Using the reweighted reconstruction of / from Q coined /, we now consider the 
following new minimization problem: 



Note that the frequency distribution of / is highly concentrated on the high 
frequencies - in other words, scaling subbands of large scales j -, and Theorem 
[Hjustifies our expectation of a very precise separation using / instead of /. Even 
more advantageous, the reduced problem (|5]) no longer involves different scales 
j, and hence can be efficiently solved by various fast numerical schemes. In our 
separation scheme, we use the same optimization method as the one used in 
MCALab to solve (|5|) now applied to a combined wavelet-shearlet dictionary. 
We refer to [7] for a detailed description of an algorithmic approach to solve ([5]); 
see also [6]. In the following subsections, we discuss the particular form of the 
matrices and ^2 which encode the wavelet and shearlet transform in the 
minimization problem (j5j) we aim to solve. 

4.4 Wavelet Transform 

Let us start with the wavelet transform. The undecimated digital wavelet trans- 
form is certainly the most fitting version of the wavelet transform for the filtering 
of data, and hence this is what we utilize also here. This transform is obtained 
by skipping the subsampling, thereby yielding an overcomplete transform, which 
in addition is shift- invariant. The redundancy factor of this transform is 3J-|- 1, 
where J is the number of decomposition levels. We refrain from further details 
and merely refer the reader to [17^ . 




4.5 Shearlet Transform 

For the shearlet transform, we employ the digital shearlet transform implemented 
by 2D convolution with discretized compactly supported shear lets, which was in- 
troduced in p!6l, see also fi3\. In the earlier work p!5], an faithful digitalization 



{W,S) = avgmin^^sW'^JWh + ll^l^^lli + M\f-W - Sg. 



(5) 
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of the continuum domain shearlet transform using compactly supported shear- 
lets generated by separable functions has been developed. However, firstly, this 
algorithmic realization allows only a limited directional selectivity due to sep- 
arability and, secondly, compactly supported shearlets generated by separable 
functions do not form a tight frame which causes an additional computational 
effort to approximate the inverse of the shearlet transform by iterative methods. 
These problems have been resolved in |16|13| by using non-separable compactly 
supported generators, and we now summarize this procedure. 

In the sequel, we will discuss the implementation strategy for computing the 
shearlet coefficients (/, V^j^/cm) only for il)j^k,m ^ ^('^^l)- The same procedure 
can be applied to shearlets in 1) except for switching the role of variables. 
Without loss of generality, let j/2 be an integer. For J > fixed, assume that 
/(^) — ^nei? /j(^)2^(/>(2^x — n), where is a 2D separable scaling function of 
the form 0i(xi)(/)i(x2) satisfying 

= ^ h{n)2(l){2x - n). (6) 

Let be a 2D separable wavelet defined by V^(xi,X2) = '0i(^i)0i(^2), where 
is a ID wavelet. Further, assume that ij) can be written as 

^l){x) = ^ w{n)2(j){2x - n), (7) 

where w{n) are 2D separable wavelet filter coefficients associated with scaling 
filter coefficients h{n). For each j > 0, define the (non-separable) shearlet gen- 
erator ij)^^^ by 

i^rio = pj-,/2{mo, (8) 

where Pi{C) = P(2^+^^i,^2) for ^ > and the trigonometric polynomial P is a 
2D fan filter (c.f. 0]). To implement {f{-),^^%\{-)) = {f{-),^r^^{S2-./2,-)), 
we make two observations: Firstly, the functions "^J^rn wavelets generated 
by refinement equations (|6]), ([7]) and (|8]). Thus, for each j, there exists an asso- 
ciated 2D wavelet filter Wj. Secondly, the shear operator 82-0 /2j^ can be faith- 
fully discretized by the digital shear operator S^_j/2j. (see |15|16j , also p!3]). 
The digital (non- separable) shearlet transform is then, using the shearlet filters 
^ik = ^2-j/2/c(^i)' defined by 

SH{fj){rm,m2) = (/j * V^^,,)(2^-^mi,2^-^/2m2) for fj e 

If downsampling by A23 is omitted, a shift invariant shearlet transform (/j * 
V^j /c)(^i' ^2) is obtained, in which case dual shearlet filters ip^ ^ can be easily 
computed by deconvolution. We then obtain the reconstruction formula 
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5 Numerical Results 

In this section, we present and discuss some numerical results of our proposed 
scheme for separating point- and curvelike features. In each experiment we com- 
pare our scheme, which is freely available in the ShearLatll toolbox, with the 
separation algorithm MCALaqj. In contrast to our algorithm, MCALab uses 
wavelets and curvelets to separate point- and curvelike components, and we re- 
fer to [7J for more details on the algorithm. 

5.1 Comparison by Visual Perception 

Aiming first at comparison by visual perception, we choose an artificial image / 
composed of two subimages P and C, where P solely contains pointlike structures 
and C different curvelike structures (see Figures [TJa) and (b)). We further add 
white Gaussian noise to / = P + C, shown in Figure [TJc). 




(a) (b) (c) 

Fig. 1. (a) P: Image of points, (b) C: Image of curves, (c) Noisy image (512 x 512). 



Our scheme consists of two parts: Preprocessing of the image as described in 
Subsection 14. 2[ followed by separation using a combined wavelet-shearlet dictio- 
nary as described in Subsections l4.3jH31 First, we will focus on the preprocessing 
step, and apply MCALab with and without our preprocessing step - due to lim- 
ited space we just mention that our scheme is similarly positively affected by 
preprocessing. Figures [2ja) and (b) then visually indicate that preprocessing 
indeed significantly improves the accuracy of separation. 

To achieve a fair comparison, we now apply both schemes to the same prepro- 
cessed image as defined in (jl]). Figures [3] (a)- (d) and Figures |ll^a)-(b) show the 
comparison results. In Figure [^c), it can be observed that the curvelet trans- 
form performs well for extracting lines due to excellent directional selectivity. 
However, some part of the curve is missed and appears in the point like part, see 
also Figure m^a). This error becomes worse with growing curvature. In contrast 
to this, compactly supported shearlets provide much better spatial localization 
than (band-Hmited) curvelets, which positively affects the capturing of local- 
ized features of the curve as illustrated in Figure [3jd) and Figure ID^b). Hence, 

^ ShearLab (Version 1.1) is available from http : //www . shearlab . org , 

^ MCALab (Version 120) is available from http : // j starck . free . f r/ j starck/Home . html [ 
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(a) (b) 

Fig. 2. (a) Pointlike image component extracted by MCALab without preprocessing, 
(b) Pointlike image component extracted by MCALab with preprocessing. 



with respect to this visual comparison, our scheme outperforms MCALab, in 
particular, when curves with large curvature are present. 



(a) 



(b) 




(d) 

Fig. 3. (a) MCALab: Pointlike component, (b) Our scheme: Pointlike component, (c) 
MCALab: Curvelike component, (d) Our scheme: Curvelike component. 



5.2 Comparison by Quantitative Measures 

We now put the comparison on more solid ground by introducing two quantita- 
tive measures for analyzing how accurate our scheme as compared to MCALab 
extracts points and curves. 

To define our first quantitative measure, let P and C be (binary) images 
containing points and curves, respectively, with image domain i? C Z^, and let 
P and C be the separated images from / = P + C + noise by the separation 
scheme to be analyzed. Letting T > 0, Bt be defined by 

Bt(/) = X{nef2:iin)>T}^ for a 2D image /, 
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and ^ be a 2D discrete Gaussian filter, we introduce the test measures 

\\9*P\\2 

and 

Using P and C as given by Figure [TJ the graphs of the error functions Mp{P) 
and MciC) for our scheme and MGALab are plotted in Figure [5] depending 
on the threshold parameter < T < 1. These figures imply that our scheme 




Fig. 5. (a) Graph of quantitative measure T ^ Mp{P){T): Our scheme (dashed curve) 
and MGALab (solid curve), (b) Graph of quantitative measure T Mp{C){T): Our 
scheme (dashed curve) and MGALab (solid curve). 



outperforms MGALab with respect to this quantitative measure. 

As our second quantitative measure, we will use the running time of each 
scheme. With respect to this comparison measure, MGALab runs 182.19 sec (30 
iterations) to produce the test results while our scheme takes 135.37 sec (15 
iterations). The running time was computed by taking the average over 10 runs. 
Again, with respect to this measure our scheme performs superior. 
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6 Application in Neurobiology 

To test the performance of our scheme on real- world images, we apply it to an 
image of a neuron generated by fluorescence microscopy (Figure [6](a)), which 
is composed of 'spines' (pointlike features) and 'dendrites' (curvelike features). 
Figures [6fb) and (c) show the extracted images containing spines and dendrites, 
respectively. 




(a) (b) (c) 

Fig. 6. (a) Image of neuron, (b) Extracted spines, (c) Extracted dendrites. 



7 Conclusion 

In this paper, we introduced a novel methodology for separating images into 
point- and curvelike features based on the new paradigm of sparse approxima- 
tion and using minimization. In contrast to other approaches, our algorithm 
utilizes a combined dictionary consisting of wavelets and shearlets, implemented 
as shift- invariant transforms, and is based on a mathematical theory. The excel- 
lent localization property of compactly supported shearlets allows shearlets to 
capture the curvehnear part very accurately and efficiently, even if the curvature 
is relatively large. Numerical results show that our scheme extracts point- and 
curvelike features more precise and uses less computing time than the state-of- 
the-art algorithm MCALab, which is based on wavelets and curvelets. 
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